#ifdef enable_superlu

#include <memory>
#include <iostream>
#include <fstream>
#include <iomanip>

#include <mpi.h>

#include "SolverSuperLU.h"

using namespace std;


SolverSuperLU::SolverSuperLU(math::SMatrix* A_, math::Vector* b_, math::Vector* x_) :
    A(A_),
    b(b_),
    x(x_)
{

    wall_time0 = new WallTime("solver init");
    wall_time1 = new WallTime("solver run");


    //wall_time0->restart();
    init();
    //wall_time0->update();


}


SolverSuperLU::~SolverSuperLU()
{
    //log_info<<"destroying solver";
    finish();
    //wall_time0->print();
    //wall_time1->print();
}


void SolverSuperLU::operator()()
{
    //wall_time1->restart();
    solve();
    //wall_time1->update();
}



void SolverSuperLU::init()
{
    nnz = A->get_nnz();

    //convert the temp "val row col" to A (compressed column format, i.e. Harwell-Boeing format)
    if ( !(a = doubleMalloc(A->get_nnz())) ) ABORT("Malloc fails for a[].");
    if ( !(asub = intMalloc(A->get_nnz())) ) ABORT("Malloc fails for asub[].");
    if ( !(xa = intMalloc(A->get_row() + 1)) ) ABORT("Malloc fails for xa[].");

    int i_ncp = 0;
    int i_nnz = 0;
    int nz_col = 0;
    int i_val = 0;


    //new algorithm
    for(int i_col = 0; i_col < A->n; i_col++)
    {
        auto& column_sm = A->data[i_col];

        nz_col = 0;
        for(auto iter = column_sm.begin(); iter != column_sm.end(); iter++)
        {
            a[i_val] = iter->second;
            asub[i_val] = iter->first;

            int i = iter->first;
            double value = iter->second;
            if(nz_col == 0)
            {
                xa[i_col] = i_val;
                nz_col    = 1;
            }
            i_val++;
        }
    }

    xa[A->get_row()] = nnz;
    //log_info<<"A and b have been structrued"<<endl;

    //================================ call superlu methods to factorize LU ==================
    /* Defaults */
    lwork = 0;
    nrhs  = 1;
    equil = YES;
    u     = 0.1;
    trans = NOTRANS;

    set_default_options(&options);
    options.Equil = equil;
    options.DiagPivotThresh = u;
    options.Trans = trans;

    options.PivotGrowth = YES;
    options.ConditionNumber = YES;

    sp_ienv(2);

    m = A->get_row();
    n = A->get_row();
    nrhs = 1;
    if ( !(rhsb = doubleMalloc(m * nrhs)) ) ABORT("Malloc fails for rhsb[].");
    if ( !(rhsx = doubleMalloc(m * nrhs)) ) ABORT("Malloc fails for rhsx[].");

    dCreate_CompCol_Matrix(&lu_A, m, n, nnz, a, asub, xa, SLU_NC, SLU_D, SLU_GE);
    dCreate_Dense_Matrix(&lu_B, m, nrhs, rhsb, m, SLU_DN, SLU_D, SLU_GE);
    dCreate_Dense_Matrix(&lu_X, m, nrhs, rhsx, m, SLU_DN, SLU_D, SLU_GE);
    xact = doubleMalloc(n * nrhs);
    ldx = n;
    dGenXtrue(n, nrhs, xact, ldx);
    dFillRHS(trans, nrhs, xact, ldx, &lu_A, &lu_B);

    if ( !(etree = intMalloc(n)) ) ABORT("Malloc fails for etree[].");
    if ( !(perm_r = intMalloc(m)) ) ABORT("Malloc fails for perm_r[].");
    if ( !(perm_c = intMalloc(n)) ) ABORT("Malloc fails for perm_c[].");
    if ( !(R = (double *) SUPERLU_MALLOC(lu_A.nrow * sizeof(double))) )
        ABORT("SUPERLU_MALLOC fails for R[].");
    if ( !(C = (double *) SUPERLU_MALLOC(lu_A.ncol * sizeof(double))) )
        ABORT("SUPERLU_MALLOC fails for C[].");
    if ( !(ferr = (double *) SUPERLU_MALLOC(nrhs * sizeof(double))) )
        ABORT("SUPERLU_MALLOC fails for ferr[].");
    if ( !(berr = (double *) SUPERLU_MALLOC(nrhs * sizeof(double))) )
        ABORT("SUPERLU_MALLOC fails for berr[].");

    /* Initialize the statistics variables. */
    StatInit(&stat);

    /* ONLY PERFORM THE LU DECOMPOSITION */
    lu_B.ncol = 0;  /* Indicate not to solve the system */

    //log_info<<"begin factorizing......"<<endl;

    dgssvx(&options, &lu_A, perm_c, perm_r, etree, equed, R, C,
           &lu_L, &lu_U, work, lwork, &lu_B, &lu_X, &rpg, &rcond, ferr, berr,
           &Glu, &mem_usage, &stat, &info);

    //log_info<<"end factorizing......"<<endl;

    //printf("LU factorization: dgssvx() returns info %d\n", info);
    //StatPrint(&stat);
    StatFree(&stat);
}


void SolverSuperLU::solve()
{
    for(int i = 0; i < b->size(); i++)
    {
        rhsb[i] = b->data[i];
    }


    /* ------------------------------------------------------------
       NOW WE SOLVE THE LINEAR SYSTEM USING THE FACTORED FORM OF A.
       ------------------------------------------------------------*/
    options.Fact = FACTORED; /* Indicate the factored form of A is supplied. */
    lu_B.ncol = nrhs;  /* Set the number of right-hand side */

    /* Initialize the statistics variables. */
    StatInit(&stat);
    dgssvx(&options, &lu_A, perm_c, perm_r, etree, equed, R, C,
           &lu_L, &lu_U, work, lwork, &lu_B, &lu_X, &rpg, &rcond, ferr, berr,
           &Glu, &mem_usage, &stat, &info);

    //printf("Triangular solve: dgssvx() returns info %d\n", info);

    for(int i = 0; i < x->size(); i++)
    {
        x->data[i] = rhsx[i];
    }


    StatFree(&stat);

}

void SolverSuperLU::finish()
{
    SUPERLU_FREE (rhsb);
    SUPERLU_FREE (rhsx);
    SUPERLU_FREE (xact);
    SUPERLU_FREE (etree);
    SUPERLU_FREE (perm_r);
    SUPERLU_FREE (perm_c);
    SUPERLU_FREE (R);
    SUPERLU_FREE (C);
    SUPERLU_FREE (ferr);
    SUPERLU_FREE (berr);
    Destroy_CompCol_Matrix(&lu_A);
    Destroy_SuperMatrix_Store(&lu_B);
    Destroy_SuperMatrix_Store(&lu_X);
    if ( lwork == 0 ) {
        Destroy_SuperNode_Matrix(&lu_L);
        Destroy_CompCol_Matrix(&lu_U);
    } else if ( lwork > 0 ) {
        SUPERLU_FREE(work);
    }
}


#endif